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ABSTRACT 



This paper compares theory and experiment for the kinetics of time-dependent 
sedimentation. We discuss non-interacting suspensions and colloids which may 
exhibit behavior similar to the one-dimensional motion of compressible gas. The 
velocity of sedimentation (or creaming) depends upon the volume fraction of the 
constituting particles and leads to Burgers-like equations for concentration profiles. 
It is shown that even the bi-dispersive system of two coupled Burgers equations has 
rich dynamics. The study of polydispersive case reveals a continuous "renormal- 
ization" of the polydispersity. We compare the Burgers system evolution with the 
experimental results on mono- and polydispersive sedimentation. The influence of 
thermal fluctuations is briefly discussed. 

PACS Numbers: 05.60, 05.40, 47.55K 
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1. Introduction 



The study of the motion of particles in a fluid goes back to Einstein and before 
that to Brown. Here we consider the effect of gravity upon the particles. If they 
are heavier than the surrounding fluid the resulting motion is called sedimenta- 
tion; if lighter it is creaming. The conventional and simplest problem involves the 
evolution of initially uniform suspension or colloid. At the bottom (of the test 
tube) there appears the sediment, and at the top - water free of particles which is 
called supernatant. Between these regions there is the original suspension itself. 
Thus, two interfaces are created, and they spread and propagate towards each 
other. The description of the bottom interface invokes high volume ratios lead- 
ing to "traffic" -like problems with possible jams and other instabilities whereby 
additional interfaces may emerge [1]. In this paper we address the evolution of 
the top interface which we can study by thinking about the dilute limit, in which 
the volume fraction of the particles is much less than 1. We apply the continuity 
equation which describes the conservation of species with concentration c{x, t) and 
current J{x,t). Such an equation is well-known in the case of very small particles 
which experience Brownian motion. It can be taken from Landau & Lifshitz [2, 
§58,59] 

dtc + div3 = (1.1) 

where the particle flux J is produced by external forces (such as gravity), and 
also depends on gradients of concentration, pressure, and temperature. In the 
experiments which we discuss below, the influence of the (other than hydrostatic) 
pressure gradient and temperature gradient is negligible, so that the particle flux 
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can be written as an expansion in the concentration gradient in the form 

J = V(c) c-D(c)Vc, (1.2) 

where the sedimentation velocity V(c) is produced by gravity and the hydrody- 
namic interactions, and the particle difFusivity D{c) is produced both by thermal 
(Brownian) motion and hydrodynamic motion. Eqs(l.l), (1.2) have not been ap- 
plied quantitatively to sedimenting systems. It is instructive to discuss the possible 
reasons. 

The kinetic coefficients V(c) and D{c) are known for small enough particles 
when the Brownian motion dominates. These mixtures are called colloids. The ve- 
locity, V(c), has been calculated in the profound paper by Batchelor [3], see also [4]. 
Batchelor also found the dependence of Brownian diffusivity upon concentration 
[5]. Batchelor's theory provides results for mono- and polydispersive suspensions. 
Yet, the task of plugging these results into Eqs(l.l), (1.2) and comparing quanti- 
tatively its evolution with real experiments has not been performed. There have 
been some qualitative results, though. It was recognized [6,7] that Eqs(l.l), (1.2) 
with c-dependent velocity resemble very closely the famous Burgers equation of 
one-dimensional compressible flow [8,2]. This equation leads naturally to shock 
waves in which the non-linear velocity dependence on concentration is balanced by 
difFusivity. We suggest here that the shock waves are the mathematical realization 
of the observed interfaces. Some non-linear aspects of this effect were studied long 
ago by Kynch [9] , who discussed the shock wave topphng in the absence of diffusiv- 
ity, and more recently by Barker and Crimson [6], and by van Saarloos and Huse 
[7]. Baker and Crimson pointed out the fact that the Burgers equation has an 
analytical solution (if the concentration dependence of difFusivity is neglected) and 
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argued its qualitative applicability to a sedimentation experiment [10]. No quan- 
titative comparison was done. The interesting paper by van Saarloos and Huse 
[7] used Burgers equation to discuss the sedimentation layers which are sometimes 
observed in the course of sedimentation. We shall discuss this phenomenon in a 
later paper [11]. 

Here we seek a comparison between the time- dependent evolution predicted 
by Burgers equation and real experimental data [12-14]. We shall also extend the 
description of sedimentation within the framework of Burgers equation to suspen- 
sions with large non-Brownian particles. Some unresolved theoretical issues arise 
in this area. For large enough particles the Brownian diffusivity is surpassed by 
the so-called hydrodynamic diffusivity. This effect has been studied by different 
experimental groups [12,14-17] but has not been fully understood theoretically, 
although a number of steps in this direction has been made [18-20]. Interesting 
numerical results have been reported by Ladd [21] for model system consisting of 
32 and 108 particles. The value of the hydrodynamic diffusivity depends upon 
unknown density-density correlation function of a suspension. This correlation 
function is not short-ranged as in the case of small enough particles (colloids) - 
the property used in the calculations of Batchelor [3,4] - and the dependence V{c) 
for suspensions is also unknown theoretically. We avoid the unresolved theoretical 
issues by making a simple order of magnitude estimate of hydrodynamic diffusivity 
which is different from previous estimates and consider the numerical prefactor as a 
fitting parameter. We also ignore the tensor-like structure of hydrodynamic diffu- 
sivity which is hopefully insignificant in the one-dimensional description offered by 
Eqs(l.l), (1.2) for concentation field depending solely upon vertical coordinate. As 
for the velocity of hindered settling, it will be desribed by an empirical dependence. 
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We would like to go further in applying the Burgers equation. This paper con- 
tains a straightforward generalization to the systems of coupled Burgers equations 
which are capable of describing realistic polydispersive suspensions. We find the 
steady-moving solutions of coupled Burgers equations and reproduce the evolu- 
tion of the suspension interface discussed originally by Smith [23] and, in implicit 
form, by Davis and Acrivos [24]. The coupled Burgers equations predict an inter- 
esting phenomenon, which we named phase shifts, it may be already observed in 
a bidispersive system. Imagine that we study the evolution of intitially homoge- 
nious bidispersive system. In a while there will be two steady-moving interfaces, 
corresponding to the two given species. Their velocities can be found by using 
Smith formulas [23] . What is the distance between these interfaces? This distance 
cannot be found by multipying the velocity difference by elapsed time. One has 
to account for the stage of "disentanglement" of the interfaces interacting via the 
non-linear coupled Burgers equations. Thus, the distance acquires the mentioned 
phase shift which is of order of combined diffusive length calculated at the disen- 
tanglement time, i.e. it is a diffusivity-dependent effect. This effect is important 
for calculation of the width of interfaces of polydispersive suspensions. On the 
other hand the quantitative definition of phase shifts is not really straightforward: 
this phenomenon becomes important when an experimentalist can spatially resolve 
the interface profile, and possible definitions contain a fair degree of uncertainlty. 
For example, one may expect a definition of the distance between the interfaces as 
a measure between the points where the concentrations of corresponding species 
diminish in, say, e times. However, in the course of sedimentation in real polydis- 
persive systems the segregation of particles improves gradually with time, and it is 
unclear which collection of species is to be traced with the above defintion. We then 
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regard the solution of the coupled Burgers equations as a necessary step in describ- 
ing the spreading of the interface. In the case of real polydispersive suspensions the 
distribution function of particle sizes a, c{a,x,t) replaces the concentration. This 
distribution function evolves in a complex way near the interface. Step by step the 
system eliminates particles with different a. We call this evolution renormalization 
since it is produced by applying many times a given rule of eliminating the fastest 
species (infinitely many times in the continuum limit). 

We introduce the Burgers equation in Section 2 and discuss the relevant kinetic 
coefficients. To show how this equation works, the experimental data of Ref.l3 for 
a colloid are compared with the analytical solution. In Section 3 we used coupled 
Burgers equations to describe polydispersive systems. We first address the unique 
steady-state motion it terms of simple analytical formulas borrowed from physics of 
one-dimensional compressible gas [2] and then simulate the bi-dispersive case nu- 
merically to investigate the evolution of the initial condition problem. Description 
of the continuous polydispersive case is presented in Section 4. In the polydisper- 
sive case sedimentation leads to the continuous renormalization of the particle size 
distribution function as one scans through the top interface. The late-time solution 
is discussed together with the effects of diffusivity. Comparison with experiments 
[12,14] is given. Section 5 contains brief discussion of influence of intrinsic thermal 
fluctuations on the sedimenting system and indicates the connection with dynam- 
ical critical phenomena. 
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2. Monodispersive sedimentation 

Depending on the particle size one distinguishes suspensions and colloids for 
a given fluid. The boundary is defined somewhat arbitrarily through a Peclet 
number, 

which describes competition between action of (say) gravity and thermal fluctu- 
ations on a spherical particle of radius a. Gravity results in particle downward 
motion with the Stokes velocity Vq, 

^ iA/3)na-^gAp ^ 2a^gAp 
dnrja 9r] ' 

due to the density difference Ap in a fluid with with viscosity r] in the presence 
of gravity acceleration g. Thermal fluctuations are the source of the Brownian 
diffusivity, 

Do=IS^. (2.3) 

where ksT is the temperature in energy units. Systems with Pe < 1 are con- 
ventionally called colloids. Thermal fluctuations are strong for colloids and lead 
to a uniform density-density correlation function which is short-ranged. In the 
opposite case the structure and kinetic properties are exclusively determined by 
the interparticle hydrodynamic interactions. For water at room temperature and 
particles less that l/im Brownian effects usually dominate, while larger particles 
form suspensions. 

Evolution of the interface which separates the region of the initial distribution 
and supernatant in colloids can be described by the concentration profile c{x, t) 



which obeys a Burgers-hke equation (see above) 

dtc = d^[V{c)c]+d^[D{c)d^c], (2.4) 

where V{c),D{c) represent the Stokes velocity and gradient diffusivity modified 
by the presence of other particles. The axis x is directed upward, opposite to the 
direction of gravity. If c(x, t) is normalized to be the volume fraction (c = ^ira^n 
with n being the number density), then 

V{c) = VoMc), , (2.5) 

where fv{c) is the hindering effect and 

fy{c) = l-kc + 0{c^), c<l. (2.6) 

According to calculations by Batchelor [3], k ~ 6.55 in the dilute limit, c <^ 1. 
Larger values of the concentration are sometimes approximated, for example, by 
the Richardson-Zaki empirical formula [22], /„(c) = (I—c/cq)'"^. (There exist other 
suggested formulas [13].) Here cq < 1 is the volume fraction at dense packing. 
The diffusivity D entering (2.4) can be written as 

Dj,{c)^DoMc), (2.7) 

with (again according to Batchelor [5], see also [25]) 

/^(c) = l-(A;-8)c+... fal + 1.45c, c < 1. (2.8) 

The nonlinear terms in the expansion (2.8) (studied experimentally [25]) lead to 
a maximum of at c ~ 0.15, so that the entire concentration dependence of the 
diffusivity is within 10% of its bare values up to the volume fractions c ~ 0.3. 
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Addressing now suspensions, we recall that the right-hand side of Eq(2.4) rep- 
resents a double expansion of the particle current in terms of small concentration 
value and small spatial gradients of concentration. In the case of large Peclet num- 
bers, we deal with strongly interacting systems, and the applicability of Eq(2.4) is 
unclear. To be more specific we do not doubt that the term V(c)c must enter the 
particle flux J in the homogenious suspension, see Eq(1.2), however the existence 
of a local (and/or fickian) expression for the inhomogeneity-induced part of flux 
has not been established. In fact, the numerical simulations by Ladd [21] indicate 
that motion of a marked particle is somewhat non-gaussian. Nevetheless, in this 
paper we use Eqs(1.2) and (2.4) to model the evolution of suspensions. Then, the 
function D{c) is conventionally called hydrodynamic diffusivity [19]. This concept 
is an attempt to get the main effect of multiparticle forces without a very detailed 
calculation and can only be justifled a posteriori, as it is sometimes the case with 
the Fokker-Planck method. 

Thus, the hydrodynamic diffusivity of suspensions, which will be denoted as 
Dfi{c), reflects the net effect of the velocity fluctuations due to interparticle hydro- 
dynamic interactions. This problem is reminiscent of those which arise in the de- 
scription of dense plasma [20]. However, hydrodynamic forces decay with particle- 
particle separation like 1/r and provide even stronger interaction as compared to 
1/r^ Coulomb forces. For Coulomb forces it is known (§27,41 in Ref.26) that plasma 
experiences collective and chaotic motions, the latter due to interparticle collisions. 
Hydrodynamic collective modes are dissipative, so the interesting question is the 
account for the collisions. Collisions in 3D plasma that are of most importance 
correspond to the small angle scattering events (§41, Ref.26, see also [27]). In 
the course of these events particles approach each other at a distance of order of 
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the Debyc screening length although the interparticle separation is much smaller 
than that. Mean "free" path is inversely proportional to the volume fraction c, 
/ ~ l/na^ a/c. Diffusivity is the product of the mean velocity fluctuations, v 
and /, -D(c) ~ av/c. This is the c-dependence which is valid for usual gases (§7,8, 
Ref.26), interacting plasma (§43, Ref.26), and it was claimed by Koch and Shaqfeh 
[20] for hydrodynamic diffusion in the dilute limit. 

Unfortunately, the pair collision approximation, which was developed by Lan- 
dau for plasma (§41, Ref.26) does not work here: for the particles with equal 
densities and of identical size there is a well-known reversibility and symmetry 
(see, for example, the concise discussion by Hinch [19]). Namely, the solution for 
viscous flow surrounding two descending spheres provides no hydrodynamical re- 
pulsion and/or attraction, and the inter-particle separation remains constant in 
time. Thus, at least triple collisions are needed, and a question immediately arises 
whether it is legitimate to use triple collisions and disregard high order ones? 
Analogous situations do arise in plasma [28]. 

Koch and Shaqfeh argued that a particle distribution without correlations is 
a solution of the Liouville continuity equations descibing all particles in a point- 
particle approximation. They continued that finite particle sizes are important to 
destroy this (say, initial) uniform distribution. Then Koch and Shaqfeh concluded 
that the most important triple collisions are those when two out of three particles 
are within 0(a) distance while the third particle is far away. We do not see any 
reason for such a conclusion. This conclusion also implies that there exists a way to 
resolve the question of higher order collisions without identifying a small parameter 
different from c <^ 1. In plasma physics this additional parameter is the condition 
of high temperature, T » e^/f (§27, Ref.26). Namely, the plasma temperature 
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must greatly exceed the Coulomb interaction at the interparticle separation f. This 
condition ensures that the Debye radius greatly exceeds f, justifies the Landau 
collision integral and so on. At present we do not think it is possible to extend 
these ideas to the physics of sedimentation. Interestingly, the numerical results by 
Ladd [21] do indicate some screening, although it is quite different from what is 
presented in [20]. 

Hinch [19] pointed an argument by Caflish that uniform distribution with the 
finite particle sizes leads to convective fiows due to density fiuctuations. Hinch 
speculated that the convection will cease at the level of mixing which is comparable 
to the interparticle separation, ac~^^^. The velocity fluctuations of the "remaining" 
convection are then of the order of the main term of the interparticle force, Vqc^/^. 
We may continue this speculation by multiplying these two values, 

Dh{c) = k'aVo, c < 1. (2.9) 

where k' is the pure number. One may expect that in the zero concentration limit 
the effect of hydrodynamic diffusivity vanishes. This does not necessarily imply 
that Dh{c) must go to zero, since the required system size and time needed to 
establish the steady-state regime both diverge as the concentration goes to zero. 
Using (2.9) one finds that the ratio Dh/Di, is of the order of the Peclet number. 
If the dependence on concentration in Eq(2.9) is incorrect, there would exist an 
additional class of mixtures which cannot be ascribed to either suspensions or 
colloids. 

Experimentalists have been analyzing hydrodynamic diffusion by two different 
techniques. One approach [12,14] is to investigate the gradient diffusion. The 
other is to trace a marked particle and evaluate the self-diffusivity [15-17]. In a 
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strongly interacting system the two approaches may, in principle, produce different 
values. Ham and Homsy [15] found that for the case of self-diffusion the factor k' 
depends upon concentration for small c < 0.01, and saturates at about k' — 5 
for larger values of c. Experiments by Nicolai et al discriminate between vertical 
and horizontal self-diffusivities and indicate /c' = 6 — 10 for the former. There is a 
possibility that the whole tensor may be inferred from experimental data by Nicolai 
et al. As for the gradient diffusion analysis it will be discussed in detail below. We 
shall estimate k' — 10 from analyzing gradient diffusion of the experiments [12,14] 
(see Section 4). This estimate is within the error bars of the experimental estimates 
of self-diffusivity. The c-dependence of was not really tested in this paper since 
both the degree of polydispersity and uncertainty of the size distribution function 
used in [12,14] are found to be too large for this purpose. 

We have introduced the Burgers equation and corresponding kinetic coefficients 
and, before we begin to use it, let us recall two approximations which are usually 
made when deriving equations (2.4)-(2.8). 

1. The particles are supposed not to interact in any fashion apart for the 
hydrodynamic interaction and excluded volume effect, i.e. these particles are hard 
spheres. Such an approximation rules out many interesting systems exhibiting 
spinodal decomposition and phase interface formation [29]. 

2. Thermal fluctuations are neglected in Eq(2.4) except for the mean-fleld 
contribution to diffusivity Di,. These are "three-dimensional" fluctuations which 
have to be averaged over the lateral plane for each height x of the sedimenting 
system. Alternatively, one can consider the 3D- version of Eq(2.4) with anisotropic 
drift in the vertical direction. As discussed by van Saarloss and Huse [7], the 
conservative noise leads only to small corrections (see also Section 5). 
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Now we attempt a detailed comparison with the experiemnt. It is well-known 
[8] that for any smooth initial distribution c{x, t) with the boundary conditions 
c{—oo,t) — Co, c{oo,t) — the late-time asymptotic is a single Burgers shock 
moving with a constant speed and width. Indeed, substituting c{x, t) = c{x + vt) 
into (2.4) and integrating from a; to oo one obtains 

D{c)^^[v-V{c)]c, (2.10) 

which for large negative x gives v — V{co). The shape of the interface can be found 
by another integration of Eq(2.10). This shape depends upon Peclet number. Using 
(2. 5)- (2. 8) we find an implicit equation for the concentration dependence on x 

(f-)Mco-c)-^.o,(-^) = ^. (.n) 

At intermediate times the experimentally resolved interface profile is time- 
dependent. We found relevant measurements in a paper by Al-Naafa and Selim 
[13] who studied the evolution of the optical density of a monodispersive colloid. 
To make the comparison, we solve Eq(2.4) with initial condition c{x, 0) = co9{x) 
to get the evolution of concentration 

cofi{x,t) 



c{x, t) 



f2(x,t) + M^,t)' 



We then try to fit the data in Fig.5, paper [13] using this formula. The result 
is shown in Fig.l. The experimental parameters are given in the figure caption. 
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Note that there are no adjustable parameters used. Initially the shape of the 
optical density profile is more narrow than the prediction of formula (2.12). After 
a transient we get the optical shape which does fit. 

We do not have any explanation of the transient. Our attempt to use more 
complex fit (see Section 4 for the description of polydispersive fit) by taking into 
account the small polydispersity of the suspension (±4.5 nm) did not result in 
significant changes. Thus the experiment both supports and denies the validity 
of Burgers equation for sedimentation. We need more experimental data to help 
resolve this issue. 

3. Bidispersive sedimentation 

We then consider the case of bimodal distribution of particle sizes. Let ci and 
C2 be the concentrations of these particles, both obeying continuity equations (1.1). 
The corresponding fiuxes Ji , J2 are 

Ji^ Vi{ci,C2) - Di{ci, 02)-^. (3.1) 

Using expansion at small concetrations we arrive to the two coupled Burgers equa- 
tions 

dtci = li9a;[(l - kuci - /ci2C2)ci] + Didlci 

dtC2 = ^25x[(l - k2lCi - k22C2)c2] + D2dlc2. (3.2) 

This system is a straightforward generalization of the mono dispersive Burgers equa- 
tion (2.4). The hindering velocity of each of the particle species depends linearly on 
both volume fractions. This dependence has been studied by Batchelor and Wen, 
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see Ref.4 where one can find method of computing the constants kn, ki2, ^21, ^22- 
The dependence of diffusivities Di, D2 upon particle concentrations will be ignored 
in this Section together with the cross-terms when a gradient of, say, c\ influences 
a flux of C2 particles. These reflnements (all vanishing with c) can wait until the 
validity of Burgers equation is understood. 

First, wc arc looking for steady-moving solutions ci^2{x — fi,2^)- Substituting 
these into equations (3.2) and integrating once we have 

= (^^1 - Vi)ci + Vxkncl - VikuciC2 (3.3) 

and an analogous second equation. Two shocks may be observed for the relevant 

boundary conditions ci^2(— 00) = coi,o2, ci^2(oo) = 0. We then have to decide which 
of the waves is faster. Without loss of generality let us assume that vi > V2 where 

ui = Vi(l - /ciicoi + /C12C02), (3.4) 

V2 = V2(l - A;2icoi + A;22Co2), (3.5) 

This leads to the first shock velocity to be vi. The substance 2 changes its con- 
centration when the first shock passes through. Let us denote this changed (we 
shall also call it renormalized) concentration as At the top of substance 1 we 
essentially have the problem of monodispersive sedimentation considered in previ- 
ous Section. Let the hindered velocity of settling be f^, as we have seen above, 
it is also the velocity of the second shock. Concentration C2 obeys the quadratic 
continuity equation in the reference frame moving with the first shock velocity vi 
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(see Ref.2, Chapter IX, and Ref.23) 



C2{V2 - vi) = 4(^2 - vi), 



(3.6) 



here V2 is 



V2^V2{1- k22cl). 



(3.7) 



Fig. 2 illustrates the steady-state geometry of the system. It can be easily shown 
that if v\ > V2 then vi > and vi — > V2 when i>i — > i>2. Consequently, no other 
steady-moving solutions exist for this system. Thus, the steady-moving solution of 
coupled Burgers system are characterized by the well-known formulas for the inter- 
face (shock) velocities [23,25,30]. It is nice to get to know that these formulas also 
describe the solutions of coupled Burgers equations. Experimental confirmation 
of shock velocities predicted by Eqs(3.6), (3.7) for bidispersive (and tridispersive) 
systems can be found in Refs.25,30. 

The coupled Burgers equations also allow to study the transient regimes prior to 
the establishing of the steady-moving shocks. Transient regimes were not reported 
experimentally, most probably due to the absence of relevant theoretical method 
of analysis, and we have no available experimental data to fit the concentration 
profiles. We hope to encourage studies of time-dependent phenomena by addressing 
physics of coupled Burgers equations. To give an example of transient behavior we 
solved the system of two coupled Burgers equations numerically. 

The two shocks when formed separate linearly in time and their asymptotical 
shapes may be found by numerical integration of Eq(3.3) with real concentration 
dependences of the diffusivities, see [5]. Both shocks have finite widths (analytical 
formulas may be obtained in a number of limits) and quickly cease to interact. 
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see Fig. 2. It must be emphasized that no matter how simple the steady-moving 
solution of the system (3.2) may seem, the transition to this solution in time may 
be rich, see Fig. 3 where a change in couphng constant k2i from 0.35 to 0.4 resulted 
in system inability to reach the dynamical steady-state within the integration time. 
Figs 2 and 3 show examples of the time dependent solutions of Eqs(3.2) which de- 
velop transient raises of concentration and even additional transient shocks (Fig.3). 
The final profiles in Fig. 2 display the two regions discuscd above: first shock at 
a; ~ 30, accompanying with renormalization of C2 from 1 up to 1.86; and the second 
shock at X f=:i 45 at the end of the integration time. Measuring rate of change of the 
separation between adjacent profiles one can verify that the shock velocity is nicely 
predicted by Eqs(3.4), (3.7). However, the separation between the shocks accounts 
for the initial condition and interaction at the stage of shock "disentanglement" as 
well as for the widths of shocks. 

It is also noteworthy that the described behavior is just an example of dynamics 
of a compressible one-dimensional gases. 

4. Poly dispersive sedimentation 

It is not hard to apply the ideas of previous Section to the real polydispersive 
case. We shall first discuss the situation where finite number N of particle sizes 
is involved and particle groups are numbered, 1 < i < A^. As before one should 
calculate the hindered velocities of sedimentation of all particle groups which are 
given by 



Vi = Vi 



j=N 

1 - J2 ^y^i] ' (4-1) 



and represent a generalization of expressions hke (3.4). Again we refer to Ref.l4 
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for the constants kij. The maximal velocity maxi(f.j) is to be found and the corre- 
sponding substance (say, io) will form the lowest (fastest) shock. Generally speak- 
ing, the motion of the edge of the size distribution is not necessarily the fastest, 
one may rather say that some substance zq will be move down first, depending 
upon the matrix kij and size (and density) distribution function. Other parti- 
cle concentrations at the top of the leader will be renormalized, according to the 
equations 



which form a system of quadratic equations. Note, that the system can be formally 
extended to i, j = io given that c^^ = 0. 

Now the entire procedure is to be repeated N times. At each level of renor- 
malization one additional substance is eliminated and others get renormalized. 
An example can be seen in Fig. 4 for the evolution of initially Gaussian-like size 
distribution function with N — 26 different particle sizes, and real parameters 
taken from paper [12]. In this particular case there is no inversion of ordering, 
i.e. larger sizes are eliminated first. The result of the calculation is a sequence of 
successively eliminated particle species and corresponding velocities. The interface 
shape becomes wider in time, although not as wide as one would obtain with- 
out renormalization. The concentration profile in the limit of vanishing diffusivity 
can be found by adding regions with renormalized concentrations at separations 
prescribed by successive velocity differences. Renormalization leads to thinner in- 
terfaces, this is the so-called phenomenon of s elf- sharpening oi the interface [12,14] 




'3 ' 




Ci{vi - Vio) 



* / * \ 



(4.3) 
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when smaller particles at the top of larger particles move faster then they would 
do if mixed with larger particles. When diffusivity is finite - a contribution to the 
interface width comes from the phase shifts. 

Before moving on to the continuous case we discuss the influence of diffusivity 
from yet a different angle. For a single Burgers shock connecting concentration 
change from c_oo to c+oo the width is given by 



where for a moment we assume the simplest possible situation: constant diffusivity 
and monodispersive case. The smaller is the concentration difference the weaker is 
the non- linear effect competing with constant diffusivity [8], and the shock width 
becomes proportionally larger. Similar phenomenon exists in the poly dispersive 
case, the larger is the number of different particle sizes which are taken into con- 
sideration - the smaller are the concentration changes, the smaller are the velocity 
differences. In addition, the widths of the successive shocks overlap. Thus, diffu- 
sivity provides a time-dependent cut-off for the achieved resolution in particle size 
step. 

This physics is reflected in the renormalization of size distribution function. 
When particles of a given size move down and are not present at some height, the 
distribution function of this size at that height becomes zero (rigorously speaking 
it is, of course, of order of exponentially small, exp{—Vx/D), diffusive corrections). 
In the limit of large N this leads to a discontinuity which moves monotonously or 
jumps depending on absence or presence of ordering inversion discussed above. Dis- 
continuous renormalization of the size distribution function implies that diffusivity 
is to be taken into account, and the existence of the diffusive time-dependent cut-off 



2D 



(4.4) 



W = 



VQk{c. 




—oo 
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regularizes the problem. At any point along the interface new structure appears as 
time goes on. The minimal resolved size step diminishes with time, and the emerg- 
ing of the fine structure of the distribution function (or segregation) continues. 
Note that certain integrated properties such as the profile of the optical density 
may be still perfectly smooth under these circumstances, it is the size distribution 
function which changes most drastically within the suspension interface. 

In the continuous case the coupled Burgers equations form an integro- 
differential equation for the evolution of the size distribution function c(a, x, t) 



dc{a, x,t) ir / \ 



dt 



oo 

— Vo{a)— c{a,x,t) — / da'k{a,a')c{a,x,t)c{a' ,x,t) + 



oo 



^J,,nM?^), (4.5, 



where kernel k{a,a') is the continuous generalization of matrix kij and diffusivity 
is a functional of c with diagonal terms of zero and first order in c and off-diagonal 
terms of first order in c. The solution of Eq(4.5) is well-defined due to the presence 
oi D. It is noteworthy that for the kernels k{a,a'),D{a,a') which dependence on 
a and a' consists of factorizable terms it may be more convenient to work with the 
Mellin transform of Eq(4.5) performed in variable a. 

Generalization of Eqs(4.2), (4.3) to the continuous case is straightforward. 
We consider the simplest case of no ordering inversion. Then at any point x, 
on the interface (at late times) there exists maximal particle size h{x) which is 
"eliminated" at this point. It is convenient to parametrize the spatial dependence 
by h{x). Then concentration becomes a function of a and 6, i.e. c(a, 6) (and 
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certainly an implicit function of x through b{x)). Particle hindered velocities are 
given by 

b 

v{a,b) ^Vo{a)-f{a,b), j{a,b) ^ 1 - J da'k{a,a')c{a'), (4.6) 



and the continuity equation reads 

c{a,b)[v{a,b) -v{b,b)] ^ c{a,b - db)[v{a,b - db) - v{b,b)], (4.7) 

where 

b-db 

v{a,b-db) = VQ{a)^l- j da'k{a,a')c{a,b - db)^. (4.8) 



Expansion of (4.6) - (4.8) to the first order to get an equation for dc{a, b)/db should 
be done carefully since c(a, b) contains a singularity at a = 6 

'^"'^^=(6^^ + ^''(«'^)' (4.9) 

where // is positive, such that integral (4.6) converges, and c^(a, b) stands for the 
regular part at a = 6. The exponent /i(6) can be related to other functions by 
using Eq(4.7), (4.8). A study not presented here shows that 

3 In 7(0,6) 

M^) = , 7^.,) ' (4-10) 

b Ba \a=b 

i.e. concentration c(a, b) becomes infinite at the point a — b. Given that the 
singularity is identified, an explicit differential equations of renormalization can be 
written for dC/db and dc^{a, b)/db. 
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Numerical integration of Eq(4.5) enables to study time-dependent evolution 
of polydispersive suspensions. For this purpose we performed simulations on a 
Sun Sparc 2 workstation to fit the experimental data by Davis and Hassen [12], 
and Lee et al [14]. The simplest possible explicit difference scheme (Euler scheme) 
already works nicely for Eq(4.5), provided that spatial and temporal steps obey 
the conditions Ax <S D/Vq, At <S min(Ax/Vb, Ax^/£)). The results cease to 
depend on numerical resolution for > 30 particle species (to model the size dis- 
tribution function) and more than 1000 spatial points for each of the species given 
that the characteristic value of the hydrodynamic diffusivity exceeds 10~^cm^/sec. 
Diffusivity of particles of size a was selected to be 



where ao is the average size. The results are presented in Figs. 4, 5 and 6. Fig. 4 
shows evolution of the size distribution function c and resembles closely the numer- 
ical results by Davis and Hassen (see their Fig. 2) if the latter are resolved in space. 
One can clearly see the appearance of singularities of the size distribution func- 
tion. Their amplitude is restricted by diffusive resolution as we discussed above. 
Fig. 5 compares simulation results with the experimental data for transmitted light 
intensity, which is given by [12] 



scaled to the experimental amplitude range /(x, 0)//(a;, oo) = 0.12. The factor 
k' was the adjustable parameter. It is interesting to note that we didn't find any 
significant dependence upon the constant k' for the values k' up to about 10. Also 
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Dh = A;'aoVb(ao), 



(4.11) 




(4.12) 



important is the observation that the interface width obtained at smaller k' (we 
tried 1 < k' < 10) is nicely comparable with the experimentally measured one. 
Smaller values of k' require higher resolution or more advanced numerical schemes 
and were not attempted. 

Our conclusion is different from the conclusion made by Davis and Hassen who 
used Eqs(4.2), (4.3) and discovered that in addition to the polydispersive width 
there is a diffusive-like contribution which they identified with the effect of gradient 
hydrodynamic diffusivity. We agree with this statement but point out that "bare" 
hydrodynamic diffusivity that enters Eqs(2.4), (4.5) doest not contribute to the 
interface width directly - it is achieved through the evolution of the equation. 
For instance, the distribution of the phase shifts influences the width. Thus, we 
think that the simulation of the coupled Burgers equation (4.5) is useful to decribe 
the results of the polydispersive sedimentation experiment. It should also be noted 
that the tails of the interface profile is sensitive to the real particle size distribution 
(which is not necessarily gaussian as assumed in Fig. 5 and in Ref.l2). 

An analogous fit was performed for the data in Fig.2 of Ref.l4. The fit is shown 
in Fig.6. Here the experimental standard deviation of particle size was only 2/3 of 
that of Ref. [12], and polydispersity alone can not account for the observed width. 
With the help of hydrodynamic diffusivity (using k' — 10) we can get curves which 
are rather similar to experimental ones. The exception is the last curve, which 
position is not described by the motion with a constant velocity. We are thankful 
to Soonchil Lee who explained that the shape and position of this curve is affected 
by the presence of the lower interface between the suspension and sediment, and 
must be ignored in our analysis [31]. 
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5. Thermal fluctuations 



Now we discuss briefly the influence of intrinsic thermal noise on sedimenta- 
tion which corresponds to Brownian motion of sedimenting particles and is only 
relevant at small Peclet numbers. One can further discriminate between two sit- 
uations. In macroscopic fluid mechanics equilibrium thermal fluctuations have 
zero correlation length. Correspondingly, different particles experience different 
Brownian drag with amplitude reflecting their sizes. One may consider adding a 
delta-correlated (in size, space and time) gaussian noise ,^(a, x, t) to the 3D gener- 
alization of Eq(4.5)(c.f. [7]). On the other side sufficiently small colloidal particles 
may be within the correlation of some thermal (especially Van-der Waals) fluctu- 
ations, then there is a correlated Brownian drag on particles of different sizes, and 
additive noise is no longer delta-correlated in particle size. This leads to complex 
behavior: for each noise-induced inhomogeneity a segregation of particles starts 
within initially uniform system. In the same manner as usual Burgers turbulence 
is an interplay between non-linearity and diffusion, we have here in addition the 
effect of size-couphng. 

Even for the monodispersive system in thermal equilibrium with average con- 
centration Co at temperature T there are fluctuations, ci = c — cq, which are 
described by Burgers equation, 

dtci = Vba^[(l - 2kcQ - kci)ci] + DAci + i. (5.1) 

Here we again neglected the ci-dependence of diffusivity, ^ is the thermal noise with 
the correlator proportional to temperature. If one averages Eq(5.1) over the lateral 
cross-section of the sedimentation column, the resulting equation is no longer 
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a closed Burgers equation but their mean-field approximations coincide. The value 
of h{x) — J'^ {ci{x,t))± is then analogous to the so-called interface height which 
was extensively studied in the framework of Burgers turbulence [32] and represents 
an example of dynamical critical phenomena. Nonequilibrium dynamics can be 
induced by changing the ambient temperature (liquid helium?) or switching off 
the interparticle interaction. 

Another question is the influence of external thermal fluctuations when ex- 
periment is performed under incomplete thermostatic conditions. This issue was 
addressed experimentally [33] and discussed theoretically [7]. The detailed study 
of this effect which results are different from Ref.7 will be published elsewhere [11]. 

6. Conclusion 

Application of (coupled) Burgers equation to describe sedimentation of (poly- 
dispersive) suspension and colloids enables one to reproduce experimental data and 
provides a basis for analyzing time-dependent sedimentation. We hope that this 
procedure will become a routine lab software and will help developing an industrial 
analysis for sedimentation processes. 
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8. Figure Captions 

Fig.l Fit of light intensity data versus tube height for a monodispersive suspension 
[13]. Curve 1 was taken after 195.2 hours after the beginning of experiment, 
curve 2 - after 374.7 hours. Experimental parameters used for the fit: a = 
65nm, cq = 0.009, Ap = 1.005 g.cm^, Vq = 0.12 x IQ-^cm/sec, T = 25^C, so 
that rj — 0.0077 g/cmsec, Dq — 4.4 x 10~*cm^/sec. Both curves are scaled 
in vertical direction to meet experimental range of Fig. 5 in Ref.l3. This is 
the experiment at small Peclet numbers, so that hydrodynamic diffusivity is 
negligible with respect to the Brownian part. 

Fig. 2 Numerical integration of two coupled Burgers equations (3.2). Parameters: 
Vi = 0.1, Va = 0.05, A;ii = ki2 = = 0.4, k2i = 0.35, L»i = L>2 = 0.01. 
Initial conditions for ci{x) and C2{x) were selected to be their individual 
Burgers shocks in the absence of coupling with small shift with respect to each 
other, c_oo = 1, c+oo — for both substances. Integration time t — 2000, 
At = 40 between profiles, ci is shown with lines, C2 with broken lines. System 
size X = 100. The sign of velocities was selected such that the shock motion 
occurs from right to left. Initial and final profiles are shown with thick lines. 
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Eqs(3.2) do not impose any contraints on amplitudes of ci,2- We can formally 
use the non-physical range ci^2 > 1- 
Fig.3 An example of complex temporal behavior. Parameters are the same as in 
Fig. 2 except for k2i = 0.4. 

Fig A Snapshot of the solution c(a, t) at time i = 36 min for the experiment 
[12]. Parameters: a = 61(±6) /xm, Ap = 1.384 g/crrt', cq = 0.05, r] — 
O.OOSSgf/cmsec. With this parameters one obtains: Vb(a) = 0.0141 cm/ sec, 
Dh{a) ^ k' X 8.5 10~^cm^/sec. The value /c' = 1 is used (see Fig.5). System 
length is 40 cm, initially gaussian size distribution with deviation 6yum is 
used. Number of species N = 26, 1600 spatial points for each particle species. 

Fig.5 Fit of raw transmitted light intensity versus time for a polydispersive sus- 
pension (Fig. 11 of Ref.l2). The intensity is measured at heights x ~ 
1.5, 6.0, 24, 32 cm. Note that we use k — 2.5; this is the experimental value 
(!) which can be easily extracted from Fig. 11 in contradiction with Table 1 
(both from Ref.l2). The other possible explanation is that Fig. 11 was, in 
fact, taken for cq = 0.02, then A; ~ 5 in agreement with Fig. 12 of Ref.l2 
and Table 1. (In both cases the conclusion made in the text remains vahd). 
The Richardson- Zaki dependence (1 — c)'^ is employed instead of 1 — kc, here 
c is the total volume fraction. The value k' — 1 is used to show that the 
hydrodynamic diffusivity is not required to account for the width of curves. 
At the top of theoretical curves one can notice small oscillations caused by 
insufficiently large number of particle species N — 26 when they get sepa- 
rated. 

Fig.6 Fit of volume fraction J da c{a,x,t) versus height x for four different times: 
332, 443, 567, and 1288 sec from left to right. Last experimental curve 



at 2611 sec is not fitted since its average displacement is not described by 
a motion with constant velocity (see text). Experimental curves are from 
Fig.2 of Ref.3. Parameters: a = 67.9(±4.0) pim, Ap = 1.57 g/cm^, rj = 
0.0164 g/cmsec, cq = 0.1. With these parameters: Vq = 0.0096 cm/ sec, 
Dfi — k'x 6.510~^ cm? /sec. The value k' — 10 is used. The Richardson-Zaki 
dependence (l — c)^ is employed with c being the total volume fraction here. 
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